Crystalline Order on a Sphere and the Generalized Thomson Problem 
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We attack generalized Thomson problems with a continuum formalism which exploits a universal 
long range interaction between defects depending on the Young modulus of the underlying lattice. 
Our predictions for the ground state energy agree with simulations of long range power law inter- 
actions of the form 1/r 7 (0 < 7 < 2) to four significant digits. The regime of grain boundaries is 
studied in the context of tilted crystalline order and the generality of our approach is illustrated 
with new results for square tilings on the sphere. 
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The Thomson problem of constructing the ground 
state of (classical) electrons interacting with a repul- 
sive Coulomb potential on a 2-spherc Jjj is almost one 
hundred years old and has many important physical 
realizations. These include multi-electron bubbles j|, 
which may be studied by capillary wave excitations, or 
the surface of liquid metal drops confined in Paul traps 
[U . Although the original Thomson problem refers to 
the ground state of spherical shells of electrons, one can 
also ask for crystalline ground states of particles interact- 
ing with other potentials. Such a generalized Thomson 
problem arises, for example, in determining the arrange- 
ments of the protein subunits which comprise the shells 
of spherical viruses [|. Here, the "particles" are clus- 
ters of protein subunits arranged on a shell. Other re- 
alizations include regular arrangements of colloidal par- 
ticles in colloidosomes proposed for encapsulation of 
active ingredients such as drugs, nutrients or living cells 
11 and fullerene patterns of carbon atoms on spheres |j] 
and other geometries p0| . An example with long range 
(logarithmic) interactions is provided by the Abrikosov 
lattice of vortices which would form at low temperatures 
in a superconducting metal shell with a large monopole 
at the center Jll[]. In practice, the "monopole" could be 
approximated by the tip of a long thin solenoid. 

Extensive numerical studies of the Thomson problem 
show that the ground state for a small number of par- 
ticles, typically M < 150, consists of twelve positive 
disclinations (the minimum number compatible with Eu- 
ler's theorem) located at the vertices of an icosahedron 
|l2| , 13 . Recent results have shown that for systems as 
small as 500 particles, however, configurations with ad- 
ditional topological defects [jl4| |l5| have lower energies 
than icosahedral ones. 

These remarkable results for the Thomson problem 
raise a number of important questions, such as the mech- 
anism behind the proliferation of defects, the nature of 
these unusual low-energy states, the universality with re- 
spect to the underlying particle potential and the gener- 
alization to more complex situations. 

A formalism suitable to address all these questions has 



been proposed recently [lq]. Disclinations are considered 
the fundamental degrees of freedom, interacting accord- 
ing to the energy [|16[ 
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where the integration is over a fixed surface with area 
element da(x) and metric gij, K is the Gaussian curva- 
ture, Y is the Young modulus in flat space and s(x) = 

^Qi${ x i Xi) is the disclination density (s(x,Xi) — 
Defects like dislocations or grain 



S(x - Xi) / y/ det(gij) 
boundaries are built from these elementary disclinations. 
Eo is the energy corresponding to a perfect defect-free 
crystal with no Gaussian curvature; Eq would be the 
ground state energy for a 2D Wigner crystal of electrons 
in the plane. Eq.(nl) restricted to a 2-sphere gives 
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where R is the radius of the sphere and x is a function of 
the geodesic distance /3y between defects with polar co- 
ordinates {6\ ; 9 j , i> j )^; x(/3) = 1 + dz Js§. 

Here 5- and 7- fold defects correspond to = +1 and —1 
respectively. In this letter we show that the continuum 
formalism embodied in Eq.(j|) implies: (a) flat space re- 
sults for elastic constants can be bootstrapped into very 
accurate quantitative calculations for generalized Thom- 
son problems, thus providing a stringent test of the va- 
lidity of this approach; (b) new results for finite length 
grain boundaries, consisting of dislocations with variable 
spacing, in the context of the 27r disclinations appear- 
ing in tilted liquid crystal phases E^]; and (c) sufficient 
power and generality to determine the ground state for 
the 8 minimal disclinations arising in square tilings of a 
sphere. 
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TABLE I: Analytical predictions (first column) for a large 
number of particles and values extrapolated from numeri- 
cal simulations (second and third columns) of the coefficient 
ai(j), as defined in Eq.(Q) for (n,n) and (n, 0) icosadeltahe- 
dral lattices. Similar accuracy holds for other values of 7. 



The Young modulus appearing in Eq.(|^) and the en- 
ergy Eq may be computed in flat space via the Ewald 
method Jl^, [lS(| . The result for M particles with long 
range pairwise interactions given by e 2 /r 7 (0 < 7 < 2) is 
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where 77, 9 and p are potential-dependent coefficients 
whose numerical values will be reported elsewhere pc[] 
and Ac is the area per particle. For M particles crys- 
tallizing on the sphere, Ac — 4irR 2 /M, and combining 
Eqs.(||) and (||) gives a large M expansion for the ground 
state energy, 



E G = 
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where 00(7) = 2 1 7 /(2 — 7) and the subleading coef- 
ficients 0^(7) depend explicitly on the potential and on 
the positions and number of disclinations. The first (non- 
extensive) term is proportional to M 2 and is usually can- 
celed by a uniform background charge for Wigner crystals 
of electrons. The coefficient a\ is a universal function of 
the positions of the defects, up to a potential-dependent 
constant. Using the results in theoretical predic- 
tions for large M for icosadeltahedral lattices [[| of type 
(n, 0) and (n, n) are given in Table [| for five values of 7. 

These predictions may now be compared with direct 
minimizations of particles on the sphere in icosadeltahe- 
dral configurations by fitting the results to Elq. (Q) . In 
Fig.|l| we plot e(M) versus 1/M for (n,0) and (n,n) 
icosadeltahedral configurations for 7 = 1.5 and 7 = 0.5, 
where 



e(M) = 
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The coefficient 01(7) is determined by the intercept in 
the M -> 00 limit (M = 10n 2 + 2 for (n, 0) lattices and 
M = 30n 2 + 2 for (n, n) lattices §). 



FIG. 1: Numerical estimate of e(M) as a function of 1/M for 
(n, 0) and (n,n) icosadeltahedral lattices with 7 = (1.5,0.5). 



(10,0) lattice 




(6,6) lattice 



FIG. 2: Strain energy distribution (red/high, blue/low) for a 
(10,0) and a (6,6) configuration. 
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FIG. 3: A 10-arm grain boundary array emerging from a +6 
defect (purple square) at the north pole. Charge +1 disclina- 
tions are red and charge —1 disclinations are yellow. 



The continuum elastic interaction between disclina- 
tions in Eq.(^) is essential to obtain the correct limiting 
behavior since it reflects contributions both from the en- 
ergy per particle in flat space as well as the energies of 
12 isolated disclinations. The defect core energy term 
in Eq.(H) contributes to the leading correction 02(7). 
We find agreement to four significant figures for (n, 0) 
icosadeltahedral lattices and to five significant figures for 
(n, n) lattices. The small residual difference in energy 
for (n, 0) and (n, n) configurations may be understood 
from the differing strain energies shown in Fig.^|. This 
small discrepancy may be attributable to a line tension 
associated with ridges of minimum or maximum strain 
connecting disclinations [ p0[ . 

We now turn to the study of grain boundaries on a 
sphere using the model described by Eq.(||). We illus- 
trate the method |L6| by considering just two 2ir defects 
(appropriate to crystals of tilted molecules) with suitable 
boundary conditions |17j] . To approximate the 2ir discli- 
nation of tilted molecules in a hemispherical crystal |l7[] 
we replace the icosahedral configuration of twelve discli- 
nations with two clusters of six 27r/6 defects at the north 
and south poles. For simplicity, we use isotropic elastic 
theory and neglect nonuniversal details near the core of 
the +27r disclination. Upon adding just one dislocation 
of Burgers vector b (i.e. a +/ — 2tt/6 disclination pair 
separated by distance b), the minimum energy in Eq.(|^) 
is achieved by a polar angle 9o(b) with b perpendicular to 
the geodesic joining the north and south poles. For small 
numbers of dislocations, the minimum energy configura- 
tion consists of two polar rings of dislocations located 
at angles 0o(b) and it — 60(b) relative to the north pole. 
The dislocations eventually organize into grain bound- 
aries centered on 9o(b), as shown in Fig.|[ Remarkably, 
no other minima were found. 

Because the global minimization just described be- 
comes computationally demanding for more than thirty 
defects, further minimizations focused on a reduced pa- 
rameter space specified by the orientations and distances 
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FIG. 4: Defect positions obtained from minimization (up tri- 
angle) and from Eq.([?j) (diamonds). 



of the grain boundaries from the two +6 defects. Fol- 
lowing B3|, the system dynamically chooses the average 
lattice spacing a = b that best accommodates the array 
structure by extremizing the energy of Eq. (||) . 

For reasonable parameter values in our continuum de- 
scription, both two antipodal +6 disclinations and the 
icosadeltahedral configurations are indeed unstable to the 
formation of grain boundaries for sufficiently large sphere 
radius[[l4| [l5| [If], 23 . In a flat monolayer, the spacing 
between dislocations in m-grain boundary arms radiating 
from a disclination of charge s is given by |21j| 
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2sin(^) 
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where b is the Burgers vector charge. The disclination 
charge is s — ^-p (p = 1 or 6, corresponding to the 
Thomson problem or tilted molecules, respectively). To 
generalize this result for symmetrical grain boundaries 
on a sphere, consider the Burgers circuit formed by an 
isosceles spherical triangle with apex angle <p g b = s/m 
at a disclination at the north pole and centered on one 
of the m-grain boundary arms p2| . If the altitude of 
this triangle is h, the net Burgers vector B defined by 
this circuit is given by the geodesic distance spanning 
the base of this triangle. A straightforward exercise in 
spherical trigonometry leads to 



cos 



cot 2 (h/R) cos 2 ((j>g b /2) + cos((/>gb) 
1 + cot 2 (h/R) cos 2 ((j) gb /2) 



(7) 



Writing B w b J dh'/l(h'), where 1(h) is a (variable) 
dislocation spacing |16| , we can invert this formula and 
thus generalize Eq.(|6|) to the sphere. 

Results comparing our minimization with Eq.(^|) are 
shown in Fig.^. Both approaches predict the same num- 
ber of dislocations within a grain and dislocation spac- 
ings which increase with 6 16 . The small discrepancies 



in the positions of the dislocations are presumably due 
to interactions between dislocations in different arms. 

The physics associated with Eq.(^) is remarkably gen- 
eral. We sketch here how the results above may be 
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adapted to a sphere tiled with a square lattice. A planar 
square lattice is described by three elastic constants (as 
opposed to the two Lame coefficients in the triangular 
case), leading to an energy 



H 



A a 



d 2 XM Q(J V , 



(8) 



where the independent elastic constants are Ai^n, An 22 
and Ai2,i2- The same derivation as that leading to Eq.(Q) 
now leads to 



H = 



1 



da(x)da(y) (k(x)-s{x)) G(x, y) (K(y)-s(y) 



(9) 

where G(x, y) = (yA 2 + 2eV?V 2 ! ) _1 , with V; the gradi- 
ent in direction i = 1 or 2, as defined by the local square 
lattice. 

tions. The elastic constants are 



The fundamental defects are now ±^ disclina- 
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The interaction energy for disclinations in a square lat- 
tice becomes equivalent to Eq.(|l]) only in the limit e = 0. 
Although details such as the critical value of R/a for the 
onset of grain boundaries will differ, we expect that the 
physics remains essentially the same. For small numbers 
of particles the ground state will consist of 8 q = +1 
disclinations. In the isotropic case (e = 0) the ground 
state is a distorted cube, with one face twisted by 45°[p0[, 
similar to tetratic liquid crystal ground states on the 
sphere Q. 

We expect similar results for geometries other than the 
sphere. For isotropic crystals on a torus with the right as- 
pect ratio, for example, one might expect 12 5-fold discli- 
nations on the outer wall (where the Gaussian curvature 
is positive) and twelve compensating 7-fold disclinations 
on the inner wall (where the Gaussian curvature is neg- 
ative) to play the role of the icosahcdral configurations 
on the sphere. As more particles are placed on the torus, 
we expect grain boundaries to emerge from these discli- 
nations. 

There are some important issues left for a future pub- 
lication p(J, such as a more detailed derivation of the 
asymptotic expansion Eq.(|]) and a reliable determina- 
tion of the optimal number of arms within a grain bound- 
ary. We hope the calculations presented here will con- 
vince the reader of the usefulness of our approach. 
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